In [1]:
using QuadGK, Plots, Roots, ForwardDiff, Interpolations, LaTeXStrings, Polynomials, PlotlyJS
using SpecialFunctions, LinearAlgebra
xb=10
RR=5
Gamma = 1.22542
masa=100
rho0 = masa/(2 *sqrt(2)*1.22542*pi*RR^3 )*0.1 ### zmienione zostalo o 0.05
G=1
gamma=0.23
delta=4*sqrt(3)*pi*0.23
plotlyjs()
#gr()
function energydensityNH(x0)
return rho0*exp(-x0^2/(2*RR^2))
end
function energydensityOS(x0)
if x0<=10
ret = rho0
elseif x0>10
ret = 0
end
return ret
end
#function Beta0(x0)
#return -asin(sqrt(8*pi*G*gamma^2*delta/x0^3 * quadgk(x->x^2*rho0*exp(-x^2/(2*RR^2)),0,x0)[1] ) )
#end
function Beta0(x0)
if x0<=10
return -asin(sqrt(8*pi*G*gamma^2*delta/x0^3 * (x0^3)/3*rho0 ) )
elseif x0>10
return -asin(sqrt(8*pi*G*gamma^2*delta/x0^3 * (10^3)/3*rho0 ) )
end
#return -asin(sqrt(8*pi*G*gamma^2*delta/x0^3 * quadgk(x->x^2*energydensityOS(x),0,x0)[1] ) )
#return -asin(sqrt(8*pi*G*gamma^2*delta/x0^3 * quadgk(x->x^2*rho0*exp(-x^2/(2*RR^2)),0,x0)[1] ) )
#return -asin(sqrt(8*pi*G*gamma^2*delta/x0^3 * rho0*(-exp(-x0^2/(2*RR^2))*RR^2*x0 +sqrt(pi/2)*RR^3*erf(x0/(sqrt(2)*RR )) ) ) )
end
function atan2(xx)
if xx<0
return atan(xx) +pi
else
return atan(xx)
end
end
function beta(lam,beta0)
return -atan2(1/(1/2*(-3*lam-2*1/tan(beta0)) ))
end
function x(lam, x0)
beta0=Beta0(x0)
return (x0*(4+9*lam^2 + 12*lam*1/tan(beta0) +4*1/(tan(beta0))^2 )^(1/3))/(2^(2/3) * (1+ 1/(tan(beta0))^2 )^(1/3) )
end
function t(lam)
return gamma*sqrt(delta)*lam
end
function lambda(t)
return t/(gamma*sqrt(delta))
end
function szift(lam, x0)
return -x(lam,x0)/(gamma*sqrt(delta))*sin(beta(lam,Beta0(x0)))*cos(beta(lam,Beta0(x0)))
end
#function xprim(lam, x0)
#beta0=Beta0(x0)
#function xx(lam)
#return (x0*(4+9*lam^2 + 12*lam*1/tan(beta0) +4*1/(tan(beta0))^2 )^(1/3))/(2^(2/3) * (1+ 1/(tan(beta0))^2 )^(1/3) )
#end
#xp = ForwardDiff.derivative(xx,lam)
#return xp
#end
function xprim(lam, x0)
return x(lam, x0)*cos(beta(lam,Beta0(x0)))*sin(beta(lam,Beta0(x0)))
end
#Plots.plot(range(-30,30,length=1000),[lam->x(lam,15),lam->xprim(lam,15)])
function tprim(lam)
return gamma*sqrt(delta)
end
Plots.plot(range(-100,100,length=1000), lam -> beta(lam, Beta0(10)) )
The WebIO Jupyter extension was not detected. See the WebIO Jupyter integration documentation for more information.
Out[1]:
In [2]:
p=Plots.plot()
for x01 in range(0,30,length=30)
Plots.plot!(range(0,100,length=500),l->x(l,x01))
end
display(p)
In [3]:
x01=10
Plots.plot(range(0,100,length=500),l->x(l,x01), color=:red)
x01=5
Plots.plot!(range(0,100,length=500),l->x(l,x01), color=:blue)
x01=11
Plots.plot!(range(0,100,length=500),l->x(l,x01))
x01=15
Plots.plot!(range(0,100,length=500),l->x(l,x01))
Out[3]:
In [4]:
function fff(lam, x0)
return -1/(1-szift(lam,x0))
end
function derfpox0(lam,x0)
function test(x0)
return fff(lam, x0)
end
return ForwardDiff.derivative(test, x0)
end
function derxpox0(lam, x0)
function test(x0)
return x(lam, x0)
end
return ForwardDiff.derivative(test, x0)
end
function derxpox0(lam, x0)
gamma1=0.23
rr=RR
x0=x0+im*0
lam=lam+im*0
rr=RR
ret=(-1 * rr ^ -1 * x0 ^ 2 * ((-2x0 + rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2)) ^ -1) ^ (2//3) * ((2x0 + -1 * rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2)) ^ -1 * (-1 * (x0 ^ 3 + 9 * G * delta * rho0 * sqrt(2) * pi ^ (3//2) * gamma1 ^ 2 * lam ^ 2 * rr ^ 3 * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1)) * exp((1 / 2) * rr ^ -2 * x0 ^ 2) + 6 * gamma1 * lam * rr * sqrt(pi) * sqrt(G * delta * rho0 * (-2x0 + rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2))) * sqrt(x0 ^ 3 + -4 * G * delta * rho0 * sqrt(2) * pi ^ (3//2) * gamma1 ^ 2 * rr ^ 3 * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) + 8 * pi * G * delta * rho0 * x0 * gamma1 ^ 2 * rr ^ 2 * exp(-(1//2) * rr ^ -2 * x0 ^ 2)) * exp(1//4 * rr ^ -2 * x0 ^ 2) + 18 * pi * G * delta * rho0 * x0 * gamma1 ^ 2 * lam ^ 2 * rr ^ 2)) ^ (-2//3) * (G * delta * rho0 * (-2x0 + rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2))) ^ (-1//2) * (x0 ^ 3 + -4 * G * delta * rho0 * sqrt(2) * pi ^ (3//2) * gamma1 ^ 2 * rr ^ 3 * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) + 8 * pi * G * delta * rho0 * x0 * gamma1 ^ 2 * rr ^ 2 * exp(-1//2 * rr ^ -2 * x0 ^ 2)) ^ (-1//2) * (-1 * rr * sqrt(G * delta * rho0 * (-2x0 + rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2))) * sqrt(x0 ^ 3 + -4 * G * delta * rho0 * sqrt(2) * pi ^ (3//2) * gamma1 ^ 2 * rr ^ 3 * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) + 8 * pi * G * delta * rho0 * x0 * gamma1 ^ 2 * rr ^ 2 * exp(-(1//2) * rr ^ -2 * x0 ^ 2)) * (6 * pi * G * delta * rho0 * gamma1 ^ 2 * lam ^ 2 + exp((1 / 2) * rr ^ -2 * x0 ^ 2)) * exp((1//4) * rr ^ -2 * x0 ^ 2) + 2 * G * delta * gamma1 * lam * rho0 * x0 * sqrt(pi) * ((x0 ^ 2 + -3 * rr ^ 2) * exp((1 / 2) * rr ^ -2 * x0 ^ 2) + 16 * pi * G * delta * rho0 * gamma1 ^ 2 * rr ^ 2) + pi * G * delta * gamma1 * lam * rho0 * sqrt(2) * rr ^ 3 * (3 * exp((1 / 2) * rr ^ -2 * x0 ^ 2) + -16 * pi * G * delta * rho0 * gamma1 ^ 2) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2)) * exp(-(5//12) * rr ^ -2 * x0 ^ 2))
return real(ret)
end
function derPRAWOpolam(lam, x0)
function test(lam)
return fff(lam,x0)*derxpox0(lam, x0)
end
return ForwardDiff.derivative(test, lam)
end
#@manipulate for lam1 in range(-40,40,length=300)
#Plots.plot(range(0.1,40,length=1000), x0-> gamma*sqrt(delta)*derfpox0(lam1,x0)+derPRAWOpolam(lam1, x0), ylims=(-10,10))
#end
Out[4]:
derPRAWOpolam (generic function with 1 method)
In [5]:
function x0odx(x0, lam)
lam1=lam
tab1=range(0.001,100,length=1000)
tab2=x.(lam1,tab1)
xx=linear_interpolation(tab2,tab1)
return xx(x0)
end
x0odx(11, 0)
#@manipulate for lam1 in range(-40,40,length=100)
##
#
#Plots.plot(range(0.1,40,length=1000), x1-> gamma*sqrt(delta)*derfpox0(lam1,x0odx(x1,lam1))
# +derPRAWOpolam(lam1, x0odx(x1,lam1)), ylims=(-10,10))
#end
Out[5]:
10.999999999999998
In [6]:
function derxpox0polam(lam, x0)
function test(lam)
return derxpox0(lam, x0)
end
return ForwardDiff.derivative(test, lam)
end
function derxpox0polam(lam, x0)
gamma1=0.23
rr=RR
x0=x0+im*0
lam=lam+im*0
rr=RR
ret=(-1 * rr ^ -1 * x0 ^ 2 * ((-2x0 + rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2)) ^ -1) ^ (2//3) * ((2x0 + -1 * rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2)) ^ -1 * (-1 * (x0 ^ 3 + 9 * G * delta * rho0 * sqrt(2) * pi ^ (3//2) * gamma1 ^ 2 * lam ^ 2 * rr ^ 3 * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1)) * exp((1 / 2) * rr ^ -2 * x0 ^ 2) + 6 * gamma1 * lam * rr * sqrt(pi) * sqrt(G * delta * rho0 * (-2x0 + rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2))) * sqrt(x0 ^ 3 + -4 * G * delta * rho0 * sqrt(2) * pi ^ (3//2) * gamma1 ^ 2 * rr ^ 3 * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) + 8 * pi * G * delta * rho0 * x0 * gamma1 ^ 2 * rr ^ 2 * exp(-1//2 * rr ^ -2 * x0 ^ 2)) * exp(1//4 * rr ^ -2 * x0 ^ 2) + 18 * pi * G * delta * rho0 * x0 * gamma1 ^ 2 * lam ^ 2 * rr ^ 2)) ^ (-2//3) * (G * delta * rho0 * (-2x0 + rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2))) ^( -1//2) * (x0 ^ 3 + -4 * G * delta * rho0 * sqrt(2) * pi ^( 3//2) * gamma1 ^ 2 * rr ^ 3 * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) + 8 * pi * G * delta * rho0 * x0 * gamma1 ^ 2 * rr ^ 2 * exp(-1//2 * rr ^ -2 * x0 ^ 2)) ^ (-1//2) * (2 * G * delta * gamma1 * rho0 * x0 * sqrt(pi) * ((x0 ^ 2 + -3 * rr ^ 2) * exp((1 / 2) * rr ^ -2 * x0 ^ 2) + 16 * pi * G * delta * rho0 * gamma1 ^ 2 * rr ^ 2) + pi * G * delta * gamma1 * rho0 * sqrt(2) * rr ^ 3 * (3 * exp((1 / 2) * rr ^ -2 * x0 ^ 2) + -16 * pi * G * delta * rho0 * gamma1 ^ 2) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2) + -12 * pi * G * delta * lam * rho0 * rr * gamma1 ^ 2 * sqrt(G * delta * rho0 * (-2x0 + rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2))) * sqrt(x0 ^ 3 + -4 * G * delta * rho0 * sqrt(2) * pi ^( 3//2) * gamma1 ^ 2 * rr ^ 3 * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) + 8 * pi * G * delta * rho0 * x0 * gamma1 ^ 2 * rr ^ 2 * exp(-1//2 * rr ^ -2 * x0 ^ 2)) * exp(1//4 * rr ^ -2 * x0 ^ 2)) * exp(-5//12 * rr ^ -2 * x0 ^ 2) + (2//3) * rr ^ -1 * x0 ^ 2 * ((-2x0 + rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2)) ^ -1) ^( 2//3) * ((2x0 + -1 * rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2)) ^ -1 * (-1 * (x0 ^ 3 + 9 * G * delta * rho0 * sqrt(2) * pi ^( 3//2) * gamma1 ^ 2 * lam ^ 2 * rr ^ 3 * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1)) * exp((1 / 2) * rr ^ -2 * x0 ^ 2) + 6 * gamma1 * lam * rr * sqrt(pi) * sqrt(G * delta * rho0 * (-2x0 + rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2))) * sqrt(x0 ^ 3 + -4 * G * delta * rho0 * sqrt(2) * pi ^ (3//2) * gamma1 ^ 2 * rr ^ 3 * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) + 8 * pi * G * delta * rho0 * x0 * gamma1 ^ 2 * rr ^ 2 * exp(-1//2 * rr ^ -2 * x0 ^ 2)) * exp(1//4 * rr ^ -2 * x0 ^ 2) + 18 * pi * G * delta * rho0 * x0 * gamma1 ^ 2 * lam ^ 2 * rr ^ 2)) ^ (-5//3) * (G * delta * rho0 * (-2x0 + rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2))) ^ (-1//2) * (2x0 + -1 * rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2)) ^ -1 * (x0 ^ 3 + -4 * G * delta * rho0 * sqrt(2) * pi ^ (3//2) * gamma1 ^ 2 * rr ^ 3 * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) + 8 * pi * G * delta * rho0 * x0 * gamma1 ^ 2 * rr ^ 2 * exp(-1//2 * rr ^ -2 * x0 ^ 2)) ^( -1//2 ) * (-1 * rr * sqrt(G * delta * rho0 * (-2x0 + rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2))) * sqrt(x0 ^ 3 + -4 * G * delta * rho0 * sqrt(2) * pi ^ (3//2) * gamma1 ^ 2 * rr ^ 3 * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) + 8 * pi * G * delta * rho0 * x0 * gamma1 ^ 2 * rr ^ 2 * exp(-1//2 * rr ^ -2 * x0 ^ 2)) * (6 * pi * G * delta * rho0 * gamma1 ^ 2 * lam ^ 2 + exp((1 / 2) * rr ^ -2 * x0 ^ 2)) * exp(1//4 * rr ^ -2 * x0 ^ 2) + 2 * G * delta * gamma1 * lam * rho0 * x0 * sqrt(pi) * ((x0 ^ 2 + -3 * rr ^ 2) * exp((1 / 2) * rr ^ -2 * x0 ^ 2) + 16 * pi * G * delta * rho0 * gamma1 ^ 2 * rr ^ 2) + pi * G * delta * gamma1 * lam * rho0 * sqrt(2) * rr ^ 3 * (3 * exp((1 / 2) * rr ^ -2 * x0 ^ 2) + -16 * pi * G * delta * rho0 * gamma1 ^ 2) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2)) * (6 * gamma1 * rr * sqrt(pi) * sqrt(G * delta * rho0 * (-2x0 + rr * sqrt(2) * sqrt(pi) * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2))) * sqrt(x0 ^ 3 + -4 * G * delta * rho0 * sqrt(2) * pi ^ (3//2) * gamma1 ^ 2 * rr ^ 3 * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) + 8 * pi * G * delta * rho0 * x0 * gamma1 ^ 2 * rr ^ 2 * exp(-1//2 * rr ^ -2 * x0 ^ 2)) * exp(1//4 * rr ^ -2 * x0 ^ 2) + 36 * pi * G * delta * lam * rho0 * x0 * gamma1 ^ 2 * rr ^ 2 + -18 * G * delta * lam * rho0 * sqrt(2) * pi ^ (3//2) * gamma1 ^ 2 * rr ^ 3 * erf((1 / 2) * x0 * sqrt(2) * rr ^ -1) * exp((1 / 2) * rr ^ -2 * x0 ^ 2)) * exp(-5//12 * rr ^ -2 * x0 ^ 2))
return real(ret)
end
Out[6]:
derxpox0polam (generic function with 1 method)
In [7]:
using DifferentialEquations
xpoczU=10
function rozwiazanie(s1,lamstar)
if s1<0
function calyproblem(du, u, p, s)
du[1] = -derxpox0(u[1], u[2])
du[2] = -gamma*sqrt(delta)
du[3] = u[3]*derxpox0polam(u[1], u[2])
end
u0 = [lamstar; xpoczU; -1/(1-szift(lamstar,xpoczU))]
sspan=(0,-s1 )
prob3 = ODEProblem(calyproblem, u0, sspan, abstol = 1e-5, reltol = 1e-5)
sol3 = solve(prob3)
return sol3(-s1)
elseif s1>=0
function calyproblem2(du, u, p, s)
du[1] = derxpox0(u[1], u[2])
du[2] = gamma*sqrt(delta)
du[3] = -u[3]*derxpox0polam(u[1], u[2])
end
u02 = [lamstar; xpoczU; -1/(1-szift(lamstar,xpoczU))]
sspan2=(0,s1+1 )
prob32 = ODEProblem(calyproblem2, u02, sspan2, abstol = 1e-5, reltol = 1e-5)
sol32 = solve(prob32)
return sol32(s1)
end
end
plotlyjs()
p=Plots.plot()
for i1 in range(-40,140,length=30)
Plots.plot!( s-> rozwiazanie(s,i1)[2], s->rozwiazanie(s,i1)[1], range(-20,10,length=500) )
end
display(p)
#Plots.plot(range(1e-15,10,length=1000),s-> sol3(s)[3])
In [8]:
function lambda1(s,lamstar)
return rozwiazanie(s,lamstar)[1]
end
function sodx0(x0)
return (x0-xpoczU)/(gamma*sqrt(delta))
end
function lamstar1(x0, lambda)
s=sodx0(x0)
tab1=range(-150,150,length=1000)
tab2=lambda1.(s,tab1)
itp=linear_interpolation(tab2,tab1)
return itp(lambda)
end
function lamstar2(x0, lambda)
s=sodx0(x0)
function test(l)
return lambda1(s,l)-lambda
end
ret=find_zero(test,-400)
return ret
end
#lamstar1(1,lambda1(-5,21))
x001=0.05
lamstar1(x001,lambda1(sodx0(x001),100))
function Ffinal(lam, x0)
if x0<xpoczU
ret = rozwiazanie(sodx0(x0),lamstar2(x0,lam))[3]
return ret
elseif x0>=xpoczU
return -1/(1-szift(lam,x0))
end
end
#Plots.plot(range(-15,35,length=500), lamstar->lambda1(sodx0(1),lamstar) )
Out[8]:
Ffinal (generic function with 1 method)
In [9]:
function derFpox0(lam,x0)
function test(x0)
return Ffinal(lam, x0)
end
return ForwardDiff.derivative(test, x0)
end
function derFpox02(lam,x0)
function test(x0)
return Ffinal(lam, x0)
end
dx0=1e-6
return (test(x0+dx0)-test(x0))/dx0
end
function derF_PRAWOpolam(lam, x0)
function test(lam)
return Ffinal(lam,x0)*derxpox0(lam, x0)
end
dlam=1e-6
return (test(lam+dlam)-test(lam))/dlam
end
#@time derFpox02(1,5)
#@time derF_PRAWOpolam(1, 1)
#lam1=20
#@time Plots.plot(range(0.1,40,length=500), x0-> gamma*sqrt(delta)*derFpox02(lam1,x0)+derF_PRAWOpolam(lam1, x0), ylims=(-10,10))
Out[9]:
derF_PRAWOpolam (generic function with 1 method)
In [10]:
function lampos(s,lamstar)
function test(s)
return rozwiazanie(s,lamstar)[1]
end
return ForwardDiff.derivative(test,s)
end
function lampos2(s,lamstar)
ds=1e-7
ret=(rozwiazanie(s+ds,lamstar)[1]-rozwiazanie(s-ds,lamstar)[1])/(2*ds)
return ret
end
@time lampos2(1,10)
0.009549 seconds (144.32 k allocations: 3.936 MiB)
Out[10]:
-0.09663934541492836
In [11]:
function Fslam(s,lamstar)
return rozwiazanie(s,lamstar)[3]
end
function FslamS(lamstar)
return -1/(1-szift(lamstar,xpoczU))
end
function intelamstarU(lamstar)
return FslamS(lamstar)*(-1)*gamma*sqrt(delta)
end
function oneoverintelamstarU(lamstar)
return 1/intelamstarU(lamstar)
end
function bisection(f, a, b, tol)
if sign(f(a)) == sign(f(b))
error("function has the same sign at given endpoints")
end
mid = (a + b)/2
while abs(f(mid)) > tol
sign(f(mid)) == sign(f(a)) ? a=mid : b=mid
mid = (a + b)/2
end
return mid
end
lamstar_zero1 = bisection(oneoverintelamstarU, -100, 0, 1e-12)
lamstar_zero2 = bisection(oneoverintelamstarU, 0, 50, 1e-12)
function calklamstarU(lamstar1)
if lamstar_zero1<lamstar1<lamstar_zero2
calk=quadgk(lamstar->intelamstarU(lamstar),0,lamstar1)[1]
elseif lamstar1<lamstar_zero1
calk=(quadgk(lamstar->intelamstarU(lamstar),0,lamstar_zero1+1e-7)[1] +
quadgk(lamstar->intelamstarU(lamstar),lamstar_zero1-1e-7,lamstar1)[1] )
elseif lamstar_zero2<lamstar1
calk=(quadgk(lamstar->intelamstarU(lamstar),0,lamstar_zero2-1e-9)[1] +
quadgk(lamstar->intelamstarU(lamstar),lamstar_zero2+1e-9,lamstar1)[1] )
end
return calk
end
function ufinalSlam(lamstar)
return 0+ calklamstarU(lamstar) #calkUslam(s,lamstar) + calklamstarU(lamstar)
end
function ufinalin(l, x0)
return 0+ calklamstarU(lamstar2(x0,l)) #calkUslam(s,lamstar) + calklamstarU(lamstar)
end
### sprawdzenie calkUslam ## JEST OK
#function dcalkUslam(s1,lamstar)
#ds1=1e-3
#return (calkUslam(s1+ds1,lamstar)-calkUslam(s1-ds1,lamstar))/(2*ds1)
#end
#star1=-10
#Plots.plot(range(-10,0,length=500), [s1->dcalkUslam(s1,star1), s1->inteUslam(s1,star1)] )
#ufinalin(1,10)
#star1=5
#Plots.plot(range(-10,0,length=500), s1->ufinalSlam(s1, star1) )
#Plots.plot(range(-120,120,length=500), lamstar-> intelamstarU(lamstar) )
Out[11]:
ufinalin (generic function with 1 method)
In [12]:
Plots.plot(range(-120,20,length=500), lamstar->ufinalSlam(lamstar))
Out[12]:
In [13]:
xpoczV=10
function rozwiazanieV(s1,lamstar)
if s1<0
function calyproblem(du, u, p, s)
du[1] = -derxpox0(u[1], u[2])
du[2] = gamma*sqrt(delta)
du[3] = u[3]*derxpox0polam(u[1], u[2])
end
u0 = [lamstar; xpoczV; 1/(1+szift(lamstar,xpoczV))]
sspan=(0,-s1 )
prob3 = ODEProblem(calyproblem, u0, sspan, abstol = 1e-5, reltol = 1e-5)
sol3 = solve(prob3)
return sol3(-s1)
elseif s1>=0
function calyproblem2(du, u, p, s)
du[1] = derxpox0(u[1], u[2])
du[2] = -gamma*sqrt(delta)
du[3] = -u[3]*derxpox0polam(u[1], u[2])
end
u02 = [lamstar; xpoczV; 1/(1+szift(lamstar,xpoczV))]
sspan2=(0,s1 )
prob32 = ODEProblem(calyproblem2, u02, sspan2, abstol = 1e-5, reltol = 1e-5)
sol32 = solve(prob32)
return sol32(s1)
end
end
plotlyjs()
p=Plots.plot()
for i1 in range(-40,140,length=30)
Plots.plot!( s-> rozwiazanieV(s,i1)[2], s->rozwiazanieV(s,i1)[1], range(-20,20,length=500) )
end
display(p)
In [14]:
function lambda1V(s,lamstar)
return rozwiazanieV(s,lamstar)[1]
end
function sodx0V(x0)
return (xpoczV-x0)/(gamma*sqrt(delta))
end
function lamstar1V(x0, lambda)
s=sodx0V(x0)
tab1=range(-550,950,length=1000)
tab2=lambda1V.(s,tab1)
itp=linear_interpolation(tab2,tab1)
return itp(lambda)
end
function lamstar2V(x0, lambda)
s=sodx0V(x0)
function test(l)
return lambda1V(s,l)-lambda
end
ret=find_zero(test,-1150)
return ret
end
#lamstar1(1,lambda1(-5,21))
x001=0.05
lamstar1V(x001,lambda1(sodx0(x001),100))
function FfinalV(lam, x0)
if x0<xpoczV
ret = rozwiazanieV(sodx0V(x0),lamstar2V(x0,lam))[3]
return ret
elseif x0>=xpoczV
return 1/(1+szift(lam,x0))
end
end
Out[14]:
FfinalV (generic function with 1 method)
In [15]:
function derFpox0V(lam,x0)
function test(x0)
return FfinalV(lam, x0)
end
return ForwardDiff.derivative(test, x0)
end
function derFpox02V(lam,x0)
function test(x0)
return FfinalV(lam, x0)
end
dx0=1e-6
return (test(x0+dx0)-test(x0))/dx0
end
function derF_PRAWOpolamV(lam, x0)
function test(lam)
return FfinalV(lam,x0)*derxpox0(lam, x0)
end
dlam=1e-6
return (test(lam+dlam)-test(lam))/dlam
end
Out[15]:
derF_PRAWOpolamV (generic function with 1 method)
In [16]:
function FslamV(s,lamstar)
return rozwiazanieV(s,lamstar)[3]
end
function FslamSV(lamstar)
return 1/(1+szift(lamstar,xpoczV))
end
function intelamstarV(lamstar)
s=0
return FslamSV(lamstar)*gamma*sqrt(delta)
end
function oneoverintelamstarV(lamstar)
return 1/intelamstarV(lamstar)
end
#Plots.plot(range(-20,30,length=500), lamstar -> oneoverintelamstarV(lamstar))
Vlamstar_zero1 = bisection(oneoverintelamstarV, 0, 7, 1e-12)
Vlamstar_zero2 = bisection(oneoverintelamstarV, 6, 105, 1e-12)
function calklamstarV(lamstar1)
if lamstar1<Vlamstar_zero1
calk=quadgk(lamstar->intelamstarV(lamstar),0,lamstar1)[1]
elseif Vlamstar_zero1<lamstar1<Vlamstar_zero2
calk=(quadgk(lamstar->intelamstarV(lamstar),0,Vlamstar_zero1-1e-9)[1] +
quadgk(lamstar->intelamstarV(lamstar),Vlamstar_zero1+1e-9,lamstar1)[1] )
elseif Vlamstar_zero2<lamstar1
calk=(quadgk(lamstar->intelamstarV(lamstar),0,Vlamstar_zero1-1e-9)[1] +
quadgk(lamstar->intelamstarV(lamstar),Vlamstar_zero1+1e-9,Vlamstar_zero2-1e-7)[1]
+quadgk(lamstar->intelamstarV(lamstar),Vlamstar_zero2+1e-7,lamstar1)[1] )
end
return calk
end
function vfinalSlam(lamstar)
return 0+ calklamstarV(lamstar) #calkUslam(s,lamstar) + calklamstarU(lamstar)
end
function vfinalin(l, x0)
return 0+ calklamstarV(lamstar2V(x0,l)) #calkUslam(s,lamstar) + calklamstarU(lamstar)
end
@time Plots.plot(range(-10,100,length=200), l->vfinalin(l, 10) )
102.284598 seconds (2.64 G allocations: 40.148 GiB, 3.79% gc time, 4.48% compilation time)
Out[16]:
In [17]:
##########################################
#####################################
######################################
######################################
function invUlam(lam,x0)
return 1/(tprim(lam)- 1/(1-szift(lam,x0))*xprim(lam, x0))
end
Plots.plot(range(-500,120,length=5000), lam->invUlam(lam,100) )
function invUlamzeros(x0)
function invUlamtest(lam)
return 1/(tprim(lam)- 1/(1-szift(lam,x0))*xprim(lam, x0))
end
return find_zeros(invUlamtest, -250,250)
end
function inteUlam(lam,x0)
return (tprim(lam)- 1/(1-szift(lam,x0))*xprim(lam, x0))
end
function calkUlam(lam,x0)
if isempty(invUlamzeros(x0))
calk = quadgk( l -> inteUlam(l, x0), 0, lam)[1]
else
zero1=invUlamzeros(x0)[1]
zero2=invUlamzeros(x0)[2]
if zero1<0 && zero2>0
if zero1<lam<zero2
calk = quadgk( l -> inteUlam(l, x0), 0, lam)[1]
elseif lam<zero1
#calk = (quadgk( l -> inteUlam(l, x0 ), 0, zero1+10e-4, rtol=10e-5 )[1]
# + quadgk( l -> inteUlam(l, x0), zero1-10e-4, lam
# , rtol=10e-5)[1] )
calk = (quadgk( l -> inteUlam(l, x0 ), 0, zero1+10e-4, )[1]
+ quadgk( l -> inteUlam(l, x0), zero1-10e-4, lam
)[1] )
#elseif zero2<lam
#calk =( quadgk( l -> inteUlam(l, x0 ), 0, zero2-10e-6, rtol=10e-7)[1]
# + quadgk( l -> inteUlam(l, x0), zero2+10e-6, lam,
# rtol=10e-7)[1] )
elseif zero2<lam
calk =( quadgk( l -> inteUlam(l, x0 ), 0, zero2-10e-7)[1]
+ quadgk( l -> inteUlam(l, x0), zero2+10e-7, lam,
)[1] )
end
elseif zero1>0 && zero2>0
if lam<zero1
calk = quadgk( l -> inteUlam(l, x0), 0, lam)[1]
elseif zero1<lam<zero2
#calk = (quadgk( l -> inteUlam(l, x0 ), 0, zero1+10e-4, rtol=10e-5 )[1]
# + quadgk( l -> inteUlam(l, x0), zero1-10e-4, lam
# , rtol=10e-5)[1] )
calk = (quadgk( l -> inteUlam(l, x0 ), 0, zero1-10e-5, )[1]
+ quadgk( l -> inteUlam(l, x0), zero1+10e-5, lam
)[1] )
#elseif zero2<lam
#calk =( quadgk( l -> inteUlam(l, x0 ), 0, zero2-10e-6, rtol=10e-7)[1]
# + quadgk( l -> inteUlam(l, x0), zero2+10e-6, lam,
# rtol=10e-7)[1] )
elseif zero2<lam
calk =(quadgk( l -> inteUlam(l, x0 ), 0, zero1-10e-5, )[1]
+ quadgk( l -> inteUlam(l, x0), zero1+10e-5, zero2-10e-5)[1] +
quadgk( l -> inteUlam(l, x0), zero2+10e-5, lam)[1] )
end
end
end
return calk
end
function xprimX0(x0)
lam1 = 0
function xxx1(x0)
lam1=0
beta0=Beta0(x0)
return (x0*(4+9*lam1^2 + 12*lam1*1/tan(beta0) +4*1/(tan(beta0))^2 )^(1/3))/(2^(2/3) * (1+ 1/(tan(beta0))^2 )^(1/3) )
end
dx0=10e-5
ret = (xxx1(x0+dx0)-xxx1(x0-dx0))/(2*dx0)
return ret
end
function xprimX000(x0)
return 1
end
function intenaX0(x0)
lam1 = 0
return -1/(1-szift(lam1, x0))*xprimX000(x0)
end
function oneoverintenaX0(x0)
return 1/intenaX0(x0)
end
function calkuodX0(x0)
if isempty(find_zeros(oneoverintenaX0, 0.1,1000))
calk = quadgk(x01 -> intenaX0(x01), 10, x0 )[1]
else
zero1=find_zeros(oneoverintenaX0, 0.1,1000)[1]
zero2=find_zeros(oneoverintenaX0, 0.1,1000)[2]
if zero1<x0<zero2
calk = quadgk(x01 -> intenaX0(x01), 10, x0 )[1]
elseif x0<zero1
calk = quadgk(x01 -> intenaX0(x01), 10, zero1+10e-8 )[1] + quadgk(x01 -> intenaX0(x01), zero1-10e-8, x0 )[1]
elseif x0>zero2
calk=quadgk(x01 -> intenaX0(x01), 10, zero2-10e-6 )[1] + quadgk(x01 -> intenaX0(x01), zero2+10e-6, x0 )[1]
end
end
return calk
end
function ufinal(lam, x0)
return calkUlam(lam,x0)+calkuodX0(x0)
end
@time ufinal(1,1)
0.559610 seconds (721.59 k allocations: 31.450 MiB, 97.05% compilation time)
Out[17]:
-1.9976469928062972
In [18]:
function invVlam(lam,x0)
return 1/(tprim(lam)+ 1/(1+szift(lam,x0))*xprim(lam, x0))
end
#Plots.plot(range(-100,900,length=5000), lam->invVlam(lam,100) )
function invVlamzeros(x0)
function invVlamtest(lam)
return 1/(tprim(lam) + 1/(1+szift(lam,x0))*xprim(lam, x0))
end
return find_zeros(invVlamtest, -300,400)
end
function inteVlam(lam,x0)
return (tprim(lam)+ 1/(1+szift(lam,x0))*xprim(lam, x0))
end
function calkVlam(lam,x0)
if isempty(invVlamzeros(x0))
calk = quadgk( l -> inteVlam(l, x0), 0, lam)[1]
else
zero1=invVlamzeros(x0)[1]
zero2=invVlamzeros(x0)[2]
if zero1<lam<zero2
calk = (quadgk( l -> inteVlam(l, x0 ), 0, zero1-10e-5)[1]
+ quadgk( l -> inteVlam(l, x0), zero1+10e-5, lam
, rtol=10e-6)[1] )
elseif lam<zero1
calk = quadgk( l -> inteVlam(l, x0), 0, lam, rtol=10e-5 )[1]
elseif zero2<lam
calk =(quadgk( l -> inteVlam(l, x0 ), 0, zero1-10e-5)[1]
+ quadgk( l -> inteVlam(l, x0), zero1+10e-5, zero2-10e-3)[1] + quadgk( l -> inteVlam(l, x0), zero2+10e-3, lam)[1] )
end
return calk
end
end
function VintenaX0(x0)
lam1 = 0
return 1/(1+szift(lam1, x0))*xprimX000(x0)
end
function VoneoverintenaX0(x0)
return 1/VintenaX0(x0)
end
#Plots.plot(range(-1,30, length=1000), x0->VintenaX0(x0))
function VcalkuodX0(x0)
calk = quadgk(x0->VintenaX0(x0), 10, x0)[1]
return calk
end
function vfinal(lam, x0)
return calkVlam(lam,x0)+VcalkuodX0(x0)
end
Out[18]:
vfinal (generic function with 1 method)
In [19]:
######## PRE COMPUTING U I V
#lamstarTAB=range(-100,100,length=5000)
#ufinalSlamTAB=ufinalSlam.(lamstarTAB)
#ufinalSlamINT=linear_interpolation(lamstarTAB,ufinalSlamTAB)
#function ufinalinINT(l, x0)
#return ufinalSlamINT(lamstar1(x0,l)) #calkUslam(s,lamstar) + calklamstarU(lamstar)
#end
####### V
#lamstarTAB=range(-100,100,length=5000)
#vfinalSlamTAB=vfinalSlam.(lamstarTAB)
#vfinalSlamINT=linear_interpolation(lamstarTAB,vfinalSlamTAB)
#function vfinalinINT(l, x0)
#return vfinalSlamINT(lamstar1(x0,l)) #calkUslam(s,lamstar) + calklamstarU(lamstar)
#end
In [20]:
kap=1/5
function atanufinal(lam, x0)
return 1/pi*atan(kap*ufinalin(lam,x0))
end
function atanufinal2(lam, x0)
return -1/pi*atan(kap*ufinalin(lam,x0))+1
end
function atanufinal3(lam, x0)
return 1/pi*atan(kap*ufinalin(lam,x0))+2
end
function atanvfinal(lam, x0)
return 1/pi*atan(kap*vfinalin(lam,x0))
end
function atanvfinal2(lam, x0)
return -1/pi*atan(kap*vfinalin(lam,x0))+1
end
function atanvfinal3(lam, x0)
return 1/pi*atan(kap*vfinalin(lam,x0))+2
end
################# TRAJEKTORIA X0=10
Ax0=10
ulamsin1=lambda1(sodx0(Ax0) ,lamstar_zero1)
ulamsin2=lambda1(sodx0(Ax0) ,lamstar_zero2)
vlamsin1=lambda1V(sodx0V(Ax0) ,Vlamstar_zero1)
vlamsin2=lambda1V(sodx0V(Ax0) ,Vlamstar_zero2)
tab1_A = range(-200,ulamsin1-10e-5,length=50)
tab2_A = range(ulamsin1+10e-3, ulamsin2-10e-7,length=50)
tab3_A = range(ulamsin2+10e-7, vlamsin1-10e-7,length=50)
tab4_A = range(vlamsin1+10e-7, vlamsin2-10e-3,length=50)
tab5_A = range(vlamsin2+10e-3, 300,length=50)
tablamA=vcat([-1000],tab1_A,tab2_A,tab3_A, tab4_A, tab5_A, [1000]);
tablamAa=tablamA
utab1_A = atanufinal.(tab1_A, Ax0)
utab2_A = atanufinal2.(tab2_A, Ax0)
utab3_A = atanufinal3.(tab3_A, Ax0)
utab4_A = atanufinal3.(tab4_A, Ax0)
utab5_A = atanufinal3.(tab5_A, Ax0)
tabutabA = vcat([-0.5],utab1_A,utab2_A, utab3_A, utab4_A, utab5_A,[2.5])
vtab1_A = atanvfinal.(tab1_A, Ax0)
vtab2_A = atanvfinal.(tab2_A, Ax0)
vtab3_A = atanvfinal.(tab3_A, Ax0)
vtab4_A = atanvfinal2.(tab4_A, Ax0)
vtab5_A = atanvfinal3.(tab5_A, Ax0)
tabvtabA = vcat([-0.5],vtab1_A,vtab2_A, vtab3_A, vtab4_A, vtab5_A,[2.5])
utrajA=linear_interpolation(tablamA, tabutabA );
vtrajA=linear_interpolation(tablamA, tabvtabA );
#### TRAJEKTORIA X0=6
Ax0=5
ulamsin1=lambda1(sodx0(Ax0) ,lamstar_zero1)
ulamsin2=lambda1(sodx0(Ax0) ,lamstar_zero2)
vlamsin1=lambda1V(sodx0V(Ax0) ,Vlamstar_zero1)
vlamsin2=lambda1V(sodx0V(Ax0) ,Vlamstar_zero2)
tab1_A = range(-200,ulamsin1-10e-1,length=50)
tab2_A = range(ulamsin1+10e-3, ulamsin2-10e-7,length=150)
tab3_A = range(ulamsin2+10e-7, vlamsin1-10e-7,length=50)
tab4_A = range(vlamsin1+10e-7, vlamsin2-10e-3,length=50)
tab5_A = range(vlamsin2+10e-3, 300,length=50)
tablamA=vcat([-1000],tab1_A,tab2_A,tab3_A, tab4_A, tab5_A, [1000]);
tablamBb=tablamA
utab1_A = atanufinal.(tab1_A, Ax0)
utab2_A = atanufinal2.(tab2_A, Ax0)
utab3_A = atanufinal3.(tab3_A, Ax0)
utab4_A = atanufinal3.(tab4_A, Ax0)
utab5_A = atanufinal3.(tab5_A, Ax0)
tabutabA = vcat([-0.5],utab1_A,utab2_A, utab3_A, utab4_A, utab5_A,[2.5])
vtab1_A = atanvfinal.(tab1_A, Ax0)
vtab2_A = atanvfinal.(tab2_A, Ax0)
vtab3_A = atanvfinal.(tab3_A, Ax0)
vtab4_A = atanvfinal2.(tab4_A, Ax0)
vtab5_A = atanvfinal3.(tab5_A, Ax0)
tabvtabA = vcat([-0.5],vtab1_A,vtab2_A, vtab3_A, vtab4_A, vtab5_A,[2.5])
utrajB=linear_interpolation(tablamA, tabutabA );
vtrajB=linear_interpolation(tablamA, tabvtabA );
################## TRAJEKTORIA X0=0.001
Ax0=0.001
ulamsin1=lambda1(sodx0(Ax0) ,lamstar_zero1)
ulamsin2=lambda1(sodx0(Ax0) ,lamstar_zero2)
vlamsin1=lambda1V(sodx0V(Ax0) ,Vlamstar_zero1)
vlamsin2=lambda1V(sodx0V(Ax0) ,Vlamstar_zero2)
tab1_A = range(-400,ulamsin1-10e-3,length=50)
tab2_A = range(ulamsin1+10e-3, ulamsin2-10e-7,length=150)
tab3_A = range(ulamsin2+10e-7, vlamsin1-10e-7,length=50)
tab4_A = range(vlamsin1+10e-7, vlamsin2-10e-3,length=50)
tab5_A = range(vlamsin2+10e-3, 300,length=50)
tablamA=vcat([-1000],tab1_A,tab2_A,tab3_A, tab4_A, tab5_A, [1000]);
tablamCc=tablamA
utab1_A = atanufinal.(tab1_A, Ax0)
utab2_A = atanufinal2.(tab2_A, Ax0)
utab3_A = atanufinal3.(tab3_A, Ax0)
utab4_A = atanufinal3.(tab4_A, Ax0)
utab5_A = atanufinal3.(tab5_A, Ax0)
tabutabA = vcat([-0.5],utab1_A,utab2_A, utab3_A, utab4_A, utab5_A,[2.5])
vtab1_A = atanvfinal.(tab1_A, Ax0)
vtab2_A = atanvfinal.(tab2_A, Ax0)
vtab3_A = atanvfinal.(tab3_A, Ax0)
vtab4_A = atanvfinal2.(tab4_A, Ax0)
vtab5_A = atanvfinal3.(tab5_A, Ax0)
tabvtabA = vcat([-0.5],vtab1_A,vtab2_A, vtab3_A, vtab4_A, vtab5_A,[2.5])
utrajC=linear_interpolation(tablamA, tabutabA );
vtrajC=linear_interpolation(tablamA, tabvtabA );
In [21]:
############# TRAJEKTORIA X0=11
## UZYTE BYLE INNE U V
function OUTatanufinal(lam, x0)
return 1/pi*atan(kap*ufinal(lam,x0))
end
function OUTatanufinal2(lam, x0)
return -1/pi*atan(kap*ufinal(lam,x0))+1
end
function OUTatanufinal3(lam, x0)
return 1/pi*atan(kap*ufinal(lam,x0))+2
end
function OUTatanvfinal(lam, x0)
return 1/pi*atan(kap*vfinal(lam,x0))
end
function OUTatanvfinal2(lam, x0)
return -1/pi*atan(kap*vfinal(lam,x0))+1
end
function OUTatanvfinal3(lam, x0)
return 1/pi*atan(kap*vfinal(lam,x0))+2
end
Ax0=11
Asin=invUlamzeros(Ax0)[1]
Asinb=invUlamzeros(Ax0)[2]
Asinc=invVlamzeros(Ax0)[1]
Asind=invVlamzeros(Ax0)[2]
tab1_A=range(-300,Asin-10e-1, length=250)
tab2_A=range(Asin+10e-1,Asinb-10e-4 , length=150)
tab3_A=range(Asinb+10e-4,Asinc-10e-4 , length=50)
tab4_A=range(Asinc+10e-4,Asind-10e-4 , length=50)
tab5_A=range(Asind+10e-4,300 , length=50)
tablamA=vcat([-1000],tab1_A,tab2_A,tab3_A, tab4_A, tab5_A, [1000]);
tablamDd=tablamA
utab1_A = OUTatanufinal.(tab1_A, Ax0)
utab2_A = OUTatanufinal2.(tab2_A, Ax0)
utab3_A = OUTatanufinal3.(tab3_A, Ax0)
utab4_A = OUTatanufinal3.(tab4_A, Ax0)
utab5_A = OUTatanufinal3.(tab5_A, Ax0)
tabutabA = vcat([-0.5],utab1_A,utab2_A, utab3_A, utab4_A, utab5_A,[2.5])
vtab1_A = OUTatanvfinal.(tab1_A, Ax0)
vtab2_A = OUTatanvfinal.(tab2_A, Ax0)
vtab3_A = OUTatanvfinal.(tab3_A, Ax0)
vtab4_A = OUTatanvfinal2.(tab4_A, Ax0)
vtab5_A = OUTatanvfinal3.(tab5_A, Ax0)
tabvtabA = vcat([-0.5],vtab1_A,vtab2_A, vtab3_A, vtab4_A, vtab5_A,[2.5])
utrajD=linear_interpolation(tablamA, tabutabA );
vtrajD=linear_interpolation(tablamA, tabvtabA );
################ TRAJEKTORIA X0=15
Ax0=15
Asin=invUlamzeros(Ax0)[1]
Asinb=invUlamzeros(Ax0)[2]
Asinc=invVlamzeros(Ax0)[1]
Asind=invVlamzeros(Ax0)[2]
tab1_A=range(-300,Asin-10e-1, length=250)
tab2_A=range(Asin+10e-1,Asinb-10e-4 , length=150)
tab3_A=range(Asinb+10e-4,Asinc-10e-4 , length=50)
tab4_A=range(Asinc+10e-4,Asind-10e-4 , length=50)
tab5_A=range(Asind+10e-4,300 , length=50)
tablamA=vcat([-1000],tab1_A,tab2_A,tab3_A, tab4_A, tab5_A, [1000]);
tablamEe=tablamA
utab1_A = OUTatanufinal.(tab1_A, Ax0)
utab2_A = OUTatanufinal2.(tab2_A, Ax0)
utab3_A = OUTatanufinal3.(tab3_A, Ax0)
utab4_A = OUTatanufinal3.(tab4_A, Ax0)
utab5_A = OUTatanufinal3.(tab5_A, Ax0)
tabutabA = vcat([-0.5],utab1_A,utab2_A, utab3_A, utab4_A, utab5_A,[2.5])
vtab1_A = OUTatanvfinal.(tab1_A, Ax0)
vtab2_A = OUTatanvfinal.(tab2_A, Ax0)
vtab3_A = OUTatanvfinal.(tab3_A, Ax0)
vtab4_A = OUTatanvfinal2.(tab4_A, Ax0)
vtab5_A = OUTatanvfinal3.(tab5_A, Ax0)
tabvtabA = vcat([-0.5],vtab1_A,vtab2_A, vtab3_A, vtab4_A, vtab5_A,[2.5])
utrajE=linear_interpolation(tablamA, tabutabA );
vtrajE=linear_interpolation(tablamA, tabvtabA );
In [22]:
#### TRAJEKTORIA X0=5
Ax0=5
ulamsin1=lambda1(sodx0(Ax0) ,lamstar_zero1)
ulamsin2=lambda1(sodx0(Ax0) ,lamstar_zero2)
vlamsin1=lambda1V(sodx0V(Ax0) ,Vlamstar_zero1)
vlamsin2=lambda1V(sodx0V(Ax0) ,Vlamstar_zero2)
tab1_A = range(-300,ulamsin1-10e-1,length=50)
tab2_A = range(ulamsin1+10e-3, ulamsin2-10e-7,length=150)
tab3_A = range(ulamsin2+10e-7, vlamsin1-10e-10,length=150)
tab4_A = vcat([vlamsin1+10e-7] ,range(vlamsin1+8, vlamsin2-10e-3,length=50))
tab5_A = range(vlamsin2+10e-3, 300,length=50)
tablamA=vcat([-1000],tab1_A,tab2_A,tab3_A, tab4_A, tab5_A, [1000]);
tablamBb=tablamA
utab1_A = atanufinal.(tab1_A, Ax0)
utab2_A = atanufinal2.(tab2_A, Ax0)
utab3_A = atanufinal3.(tab3_A, Ax0)
utab4_A = atanufinal3.(tab4_A, Ax0)
utab5_A = atanufinal3.(tab5_A, Ax0)
tabutabA = vcat([-0.5],utab1_A,utab2_A, utab3_A, utab4_A, utab5_A,[2.5])
vtab1_A = atanvfinal.(tab1_A, Ax0)
vtab2_A = atanvfinal.(tab2_A, Ax0)
vtab3_A = atanvfinal.(tab3_A, Ax0)
vtab4_A = atanvfinal2.(tab4_A, Ax0)
vtab5_A = atanvfinal3.(tab5_A, Ax0)
tabvtabA = vcat([-0.5],vtab1_A,vtab2_A, vtab3_A, vtab4_A, vtab5_A,[2.5])
utrajB=linear_interpolation(tablamA, tabutabA );
vtrajB=linear_interpolation(tablamA, tabvtabA );
In [23]:
[ulamsin1, ulamsin2, vlamsin1, vlamsin2]
Out[23]:
4-element Vector{Float64}:
-110.63904313261818
-5.236859665981674
6.757238501613085
116.048470561546
In [24]:
[ulamsin1, ulamsin2, vlamsin1, vlamsin2]
Out[24]:
4-element Vector{Float64}:
-110.63904313261818
-5.236859665981674
6.757238501613085
116.048470561546
In [25]:
function lamsin1(x0)
ulamsin1=lambda1(sodx0(x0) ,lamstar_zero1)
return ulamsin1
end
function lamsin2(x0)
ulamsin2=lambda1(sodx0(x0) ,lamstar_zero2)
return ulamsin2
end
function lamsin3(x0)
vlamsin1=lambda1V(sodx0V(x0) ,Vlamstar_zero1)
return vlamsin1
end
function lamsin4(x0)
vlamsin2=lambda1V(sodx0V(x0) ,Vlamstar_zero2)
return vlamsin2
end
Out[25]:
lamsin4 (generic function with 1 method)
In [26]:
#pv
function minusNzeros(x0)
function test(lam)
return 1-szift(lam, x0)
end
lamtab1=range(-1000,1750,length=40000)
appalamtab=[]
for i1 in lamtab1
push!(appalamtab, test(i1))
end
miejsce=lamtab1[findmin(appalamtab)[2]]
return find_zeros(test, (miejsce-150,miejsce+150))
end
function plusNzeros(x0)
function test(lam)
return 1+szift(lam, x0)
end
lamtab1=range(-1000,1750,length=40000)
appalamtab=[]
for i1 in lamtab1
push!(appalamtab, test(i1))
end
miejsce=lamtab1[findmin(appalamtab)[2]]
return find_zeros(test, (miejsce-150,miejsce+450))
end
st=2.97
Plots.plot(range(0.1,25,length=500), x0->lamsin1(x0), linewidth=2)
Plots.plot!(range(0.1,25,length=500), x0->lamsin2(x0), linewidth=2)
Plots.plot!(range(0.1,25,length=500), x0->lamsin3(x0), linewidth=2)
Plots.plot!(range(0.1,25,length=500), x0->lamsin4(x0), linewidth=2)
Plots.plot!(range(st,25,length=500), x0->minusNzeros(x0)[1], color=:black)
Plots.plot!(range(st,25,length=500), x0->minusNzeros(x0)[2], color=:black)
Plots.plot!(range(st,25,length=500), x0->plusNzeros(x0)[1], color=:black)
Plots.plot!(range(st,25,length=500), x0->plusNzeros(x0)[2], color=:black, xlims=(0,20))
Out[26]:
In [27]:
function tescik(x0)
return minusNzeros(x0)[1] - lamsin2(x0)
end
#tescik(5.25)
przec1 = find_zero(tescik, (st,10) )
Out[27]:
5.061379984358704
In [28]:
#### APPARENT HORIZONS
function app_v_dol1(x0)
return atanvfinal(minusNzeros(x0)[1],x0)
end
function app_v_dol2(x0)
return atanvfinal(minusNzeros(x0)[2],x0)
end
function app_u_dol1(x0)
return atanufinal2(minusNzeros(x0)[1],x0)
end
function app_u_dol2(x0)
return atanufinal2(minusNzeros(x0)[2],x0)
end
function app_u_dol1W(x0)
return atanufinal3(minusNzeros(x0)[1],x0)
end
function app_u_dol2W(x0)
return atanufinal3(minusNzeros(x0)[2],x0)
end
### APP DOL INSIDE
x0ran=range(przec1,9.999,length=150)
app_vtab1=app_v_dol1.(x0ran)
app_utab1=app_u_dol1.(x0ran)
linapp_v = linear_interpolation(x0ran, app_vtab1)
linapp_u = linear_interpolation(x0ran, app_utab1)
Plots.plot( x0->linapp_v(x0) - linapp_u(x0), x0->linapp_v(x0) + linapp_u(x0), x0ran, color=:magenta )
x0ran2 = range(st, przec1, length=150)
app_vtab1_przec = app_v_dol1.(x0ran2)
app_utab1_przec = app_u_dol1W.(x0ran2)
linapp_v_przec=linear_interpolation(x0ran2,app_vtab1_przec)
linapp_u_przec=linear_interpolation(x0ran2, app_utab1_przec)
Plots.plot!( x0->linapp_v_przec(x0) - linapp_u_przec(x0), x0->linapp_v_przec(x0) + linapp_u_przec(x0), x0ran2, color=:magenta )
Plots.plot!( [linapp_v(przec1) - linapp_u(przec1),linapp_v_przec(przec1) - linapp_u_przec(przec1) ],
[linapp_v(przec1) + linapp_u(przec1),linapp_v_przec(przec1) + linapp_u_przec(przec1) ], color=:magenta )
x0ran3= range(st, 10, length=150)
app_vtab2=app_v_dol2.(x0ran3)
app_utab2=app_u_dol2W.(x0ran3)
linapp_v2 = linear_interpolation(x0ran3, app_vtab2)
linapp_u2 = linear_interpolation(x0ran3, app_utab2)
Plots.plot!( x0->linapp_v2(x0) - linapp_u2(x0), x0->linapp_v2(x0) + linapp_u2(x0), x0ran3, color=:magenta )
Plots.plot!( [linapp_v_przec(st) - linapp_u_przec(st), linapp_v2(st) - linapp_u2(st)],
[linapp_v_przec(st) + linapp_u_przec(st), linapp_v2(st) + linapp_u2(st)], color=:magenta )
Plots.plot!( [OUTatanvfinal(lamsin1(10) ,10) - 0.5 , 0.5 -0.5 ], [OUTatanvfinal(lamsin1(10) ,10) + 0.5 , 0.5 +0.5 ], color=:magenta )
### APP DOL OUTSIDE
#x0ran=range(11,50,length=150)
#app_vtab1=app_v_dol1.(x0ran)
#app_utab1=app_u_dol1.(x0ran)
#linapp_vOUT = linear_interpolation(x0ran, app_vtab1)
#linapp_uOUT = linear_interpolation(x0ran, app_utab1)
#Plots.plot!( x0->linapp_v(x0) - linapp_u(x0), x0->linapp_v(x0) + linapp_u(x0), x0ran, color=:blue )
Out[28]:
In [29]:
### APP GORNY
function tescik2(x0)
return plusNzeros(x0)[2] - lamsin3(x0)
end
#tescik(5.25)
przec2 = find_zero(tescik2, (st,10) )
function upAPPv(x0)
return atanvfinal(plusNzeros(x0)[1],x0)
end
function uppAPPu(x0)
return atanufinal3(plusNzeros(x0)[1],x0)
end
function upAPPvB(x0)
return atanvfinal(plusNzeros(x0)[2],x0)
end
function uppAPPuB(x0)
return atanufinal3(plusNzeros(x0)[2],x0)
end
function upAPPvBW(x0)
return atanvfinal2(plusNzeros(x0)[2],x0)
end
function uppAPPuBW(x0)
return atanufinal3(plusNzeros(x0)[2],x0)
end
up_x0ran=range(st,9.999, length=150)
uptabav1= upAPPv.(up_x0ran)
uptabau1= uppAPPu.(up_x0ran)
linupav1 = linear_interpolation(up_x0ran, uptabav1)
linupau1 = linear_interpolation(up_x0ran, uptabau1)
## APP wyzej
up_x0ran2=range(st, przec2, length=150)
uptabavB= upAPPvB.(up_x0ran2)
uptabauB= uppAPPuB.(up_x0ran2)
linupavB = linear_interpolation(up_x0ran2, uptabavB)
linupauB = linear_interpolation(up_x0ran2, uptabauB)
## APP PO PRZECIECIU
up_x0ran3=vcat(range(przec2,przec2+1e-3 , length=75), range(4.1,9.999,length=75))
uptabavBW= upAPPvBW.(up_x0ran3)
uptabauBW= uppAPPuBW.(up_x0ran3)
linupavBW = linear_interpolation(up_x0ran3, uptabavBW)
linupauBW = linear_interpolation(up_x0ran3, uptabauBW)
Plots.plot( x0-> linupav1(x0) -linupau1(x0), x0-> linupav1(x0) +linupau1(x0), up_x0ran, color=:darkviolet )
Plots.plot!( x0-> linupavB(x0) -linupauB(x0), x0-> linupavB(x0) +linupauB(x0), up_x0ran2, color=:darkviolet )
Plots.plot!( x0-> linupavBW(x0) -linupauBW(x0), x0-> linupavBW(x0) +linupauBW(x0), up_x0ran3, color=:darkviolet )
Plots.plot!( [linupav1(st)-linupau1(st) ,linupavB(st)-linupauB(st) ], [linupav1(st)+linupau1(st) ,linupavB(st)+linupauB(st) ], color=:darkviolet )
Plots.plot!( [linupav1(st)-linupau1(st) , linupavB(st)-linupauB(st) ], [linupav1(st)+linupau1(st) ,linupavB(st)+linupauB(st) ],
color=:darkviolet )
Plots.plot!( [linupavBW(przec2)-linupauBW(przec2) , linupavB(przec2)-linupauB(przec2) ],
[linupavBW(przec2)+linupauBW(przec2) ,linupavB(przec2)+linupauB(przec2) ],
color=:darkviolet )
Out[29]:
In [30]:
Plots.plot( x0-> linupav1(x0) -linupau1(x0), x0-> linupav1(x0) +linupau1(x0), up_x0ran, color=:darkviolet )
Plots.plot!( x0-> linupavB(x0) -linupauB(x0), x0-> linupavB(x0) +linupauB(x0), up_x0ran2, color=:darkviolet )
Plots.plot!( x0-> linupavBW(x0) -linupauBW(x0), x0-> linupavBW(x0) +linupauBW(x0), up_x0ran3, color=:darkviolet )
Plots.plot!( [linupav1(st)-linupau1(st) , linupavB(st)-linupauB(st) ], [linupav1(st)+linupau1(st) ,linupavB(st)+linupauB(st) ],
color=:darkviolet )
Plots.plot!( [linupavBW(przec2)-linupauBW(przec2) , linupavB(przec2)-linupauB(przec2) ],
[linupavBW(przec2)+linupauBW(przec2) ,linupavB(przec2)+linupauB(przec2) ],
color=:darkviolet )
Out[30]:
In [31]:
## OBA APPARENT
Plots.plot( x0->linapp_v(x0) - linapp_u(x0), x0->linapp_v(x0) + linapp_u(x0), x0ran, color=:magenta )
Plots.plot!( x0->linapp_v_przec(x0) - linapp_u_przec(x0), x0->linapp_v_przec(x0) + linapp_u_przec(x0), x0ran2, color=:magenta )
Plots.plot!( x0->linapp_v2(x0) - linapp_u2(x0), x0->linapp_v2(x0) + linapp_u2(x0), x0ran3, color=:magenta )
Plots.plot!( [linapp_v(przec1) - linapp_u(przec1),linapp_v_przec(przec1) - linapp_u_przec(przec1) ],
[linapp_v(przec1) + linapp_u(przec1),linapp_v_przec(przec1) + linapp_u_przec(przec1) ], color=:magenta )
Plots.plot!( [linapp_v_przec(st) - linapp_u_przec(st), linapp_v2(st) - linapp_u2(st)],
[linapp_v_przec(st) + linapp_u_przec(st), linapp_v2(st) + linapp_u2(st)], color=:magenta )
Plots.plot!( [OUTatanvfinal(lamsin1(10) ,10) - 0.5 , 0.5 -0.5 ], [OUTatanvfinal(lamsin1(10) ,10) + 0.5 , 0.5 +0.5 ], color=:magenta )
Plots.plot!( [OUTatanvfinal(lamsin2(10) ,10) - 1.5 , 0.5 -1.5 ], [OUTatanvfinal(lamsin2(10) ,10) + 1.5 , 0.5 +1.5 ], color=:magenta )
Plots.plot!([OUTatanvfinal(lamsin2(10) ,10) - 1.5 , linapp_v2(9.999) - linapp_u2(9.999) ],
[OUTatanvfinal(lamsin2(10) ,10) + 1.5 , linapp_v2(9.999) + linapp_u2(9.999) ], color=:magenta)
Plots.plot!( [OUTatanvfinal(lamsin1(10) ,10) - 0.5 , linapp_v(9.999) - linapp_u(9.999) ],
[OUTatanvfinal(lamsin1(10) ,10) + 0.5 , linapp_v(9.999) + linapp_u(9.999) ], color=:magenta )
## DOLNY
Plots.plot!( x0-> linupav1(x0) -linupau1(x0), x0-> linupav1(x0) +linupau1(x0), up_x0ran, color=:darkviolet )
Plots.plot!( x0-> linupavB(x0) -linupauB(x0), x0-> linupavB(x0) +linupauB(x0), up_x0ran2, color=:darkviolet )
Plots.plot!( x0-> linupavBW(x0) -linupauBW(x0), x0-> linupavBW(x0) +linupauBW(x0), up_x0ran3, color=:darkviolet )
Plots.plot!( [linupav1(st)-linupau1(st) , linupavB(st)-linupauB(st) ], [linupav1(st)+linupau1(st) ,linupavB(st)+linupauB(st) ],
color=:darkviolet )
Plots.plot!( [linupavBW(przec2)-linupauBW(przec2) , linupavB(przec2)-linupauB(przec2) ],
[linupavBW(przec2)+linupauBW(przec2) ,linupavB(przec2)+linupauB(przec2) ],
color=:darkviolet )
Out[31]:
In [32]:
## TIMELIKE SINGULARITY
function trajx(T)
return x(T/(gamma*sqrt(delta)),10)
end
M= 4/3*pi*xb^3*rho0
ff = x -> 1 - 2*G*M/x + 4*G^2*M^2*gamma^2*delta/x^4
xs1=find_zeros(ff,(0,150) )[1]
xs2=find_zeros(ff,(0,150) )[2]
inte =x-> 1/(ff(x) )
xszift = x-> sqrt(2*G*M*(x^3-2*G*M*gamma^2*delta)/x^4)
function calkX(x)
if x<xs1
ret = quadgk( x1->inte(x1), 1e-15, x)[1]
elseif xs1<x<xs2
ret = quadgk( x1->inte(x1), 0, xs1-1e-7)[1] + quadgk( x1->inte(x1), xs1+1e-7, x)[1]
elseif xs2<x
ret = quadgk( x1->inte(x1), 0, xs1-1e-7)[1] + quadgk( x1->inte(x1), xs1+1e-12, xs2-1e-7)[1] + quadgk( x1->inte(x1), xs2+1e-7, x)[1]
end
return ret
end
intet = T-> 1/ff(trajx(T))
tsin1 = find_zeros(intet, (-100,100) )[1]
tsin2 = find_zeros(intet, (-100,100) )[2]
tsin3 = find_zeros(intet, (-100,100) )[3]
tsin4 = find_zeros(intet, (-100,100) )[4]
function calkt(T)
if T<tsin1
ret= quadgk( T1->intet(T1), 0, tsin1+1e-7 )[1] + quadgk( T1->intet(T1), tsin1-1e-7, T )[1]
elseif tsin1<T<tsin2
ret = quadgk( T1->intet(T1), 0, T )[1]
elseif tsin2<T<tsin3
ret= quadgk( T1->intet(T1), 0, tsin2-1e-7 )[1] + quadgk( T1->intet(T1), tsin2+1e-7, T)[1]
elseif tsin3<T<tsin4
ret = (quadgk( T1->intet(T1), 0, tsin2-1e-7 )[1] + quadgk( T1->intet(T1), tsin2+1e-7, tsin3-1e-7 )[1]
+quadgk( T1->intet(T1), tsin3+1e-7, T )[1] )
elseif tsin4<T
ret = (quadgk( T1->intet(T1), 0, tsin2-1e-7 )[1] + quadgk( T1->intet(T1), tsin2+1e-7, tsin3-1e-7 )[1]+
quadgk( T1->intet(T1), tsin3+1e-7, tsin4-1e-7 )[1]
+quadgk( T1->intet(T1), tsin4+1e-7, T)[1] )
end
return ret
end
function usin(t,x)
return -1/pi*atan( calkt(t) - calkX(x) ) +1
end
function vsin(t,x)
return -1/pi*atan( calkt(t) + calkX(x) ) +1
end
x1=0.5
Plots.plot( t-> vsin(t,x1)-usin(t,x1), t-> vsin(t,x1)+usin(t,x1), range(-10,10,length=100), color=:black )
Out[32]:
In [33]:
######### POWIERZCHNIA POCZATKOWA
in_przec1=find_zeros(lamsin2,(0,10))[1]
in_przec2=61.5558
#Plots.plot(x0->atanvfinal(0,x0)-atanufinal3(0,x0), x0->atanvfinal(0,x0)+atanufinal3(0,x0), range(1e-5,in_przec1,length=100), color=:black)
#Plots.plot!(x0->atanvfinal(0,x0)-atanufinal2(0,x0), x0->atanvfinal(0,x0)+atanufinal2(0,x0), range(in_przec1,10,length=100), color=:black)
#Plots.plot!(x0->OUTatanvfinal(0,x0)-OUTatanufinal2(0,x0), x0->OUTatanvfinal(0,x0)+OUTatanufinal2(0,x0), range(10,in_przec2-1e-1,length=100),
# color=:black)
#Plots.plot!(x0->OUTatanvfinal(0,x0)-OUTatanufinal(0,x0), x0->OUTatanvfinal(0,x0)+OUTatanufinal(0,x0),
#range(in_przec2+1e-1,120,length=100), color=:black)
intabx01=range(1e-5,in_przec1,length=100)
intabx02=range(in_przec1,10-1e-15,length=100)
intabx03=range(10+1e-15,in_przec2-1e-1,length=100)
intabx04=range(in_przec2+1e-1,120,length=100)
intabx0=vcat(intabx01,intabx02,intabx03,intabx04, [1000])
intabv1=atanvfinal.(0,intabx01)
intabv2=atanvfinal.(0,intabx02)
intabv3=OUTatanvfinal.(0,intabx03)
intabv4=OUTatanvfinal.(0,intabx04)
intabu1=atanufinal3.(0,intabx01)
intabu2=atanufinal2.(0,intabx02)
intabu3=OUTatanufinal2.(0,intabx03)
intabu4=OUTatanufinal.(0,intabx04)
intabv=vcat(intabv1,intabv2,intabv3,intabv4, [0.5])
intabu=vcat(intabu1,intabu2,intabu3,intabu4, [-0.5])
innv=linear_interpolation( intabx0, intabv)
innu=linear_interpolation( intabx0, intabu)
Plots.plot( x0-> innv(x0)-innu(x0), x0-> innv(x0)+innu(x0), intabx0, color=:black)
┌ Warning: Duplicated knots were deduplicated. Use Interpolations.deduplicate_knots!(knots) explicitly to avoid this warning. │ k1 = │ 401-element Vector{Float64}: │ 1.0e-5 │ 0.07401944551836069 │ 0.14802889103672137 │ 0.22203833655508207 │ 0.29604778207344273 │ ⋮ │ 118.23199393939394 │ 118.82132929292929 │ 119.41066464646465 │ 120.0 │ 1000.0 └ @ Interpolations C:\Users\mihalbobula\.julia\packages\Interpolations\91PhN\src\gridded\gridded.jl:77
Out[33]:
In [37]:
############ DIAGRAM
#plotlyjs()
gr()
######## APPPARENT DOLNY
lwa=1.4
cap=:magenta
Plots.plot( x0->linapp_v(x0) - linapp_u(x0), x0->linapp_v(x0) + linapp_u(x0), x0ran, color=cap,
label=L"\mathrm{first} \, \mathrm{AH}", lw=lwa )
Plots.plot!( x0->linapp_v_przec(x0) - linapp_u_przec(x0), x0->linapp_v_przec(x0) + linapp_u_przec(x0), x0ran2, color=cap, label=false, lw=lwa )
Plots.plot!( x0->linapp_v2(x0) - linapp_u2(x0), x0->linapp_v2(x0) + linapp_u2(x0), x0ran3, color=cap, label=false, lw=lwa )
Plots.plot!( [linapp_v(przec1) - linapp_u(przec1),linapp_v_przec(przec1) - linapp_u_przec(przec1) ],
[linapp_v(przec1) + linapp_u(przec1),linapp_v_przec(przec1) + linapp_u_przec(przec1) ], color=cap, label=false, lw=lwa )
Plots.plot!( [linapp_v_przec(st) - linapp_u_przec(st), linapp_v2(st) - linapp_u2(st)],
[linapp_v_przec(st) + linapp_u_przec(st), linapp_v2(st) + linapp_u2(st)], color=cap, label=false, lw=lwa )
Plots.plot!( [OUTatanvfinal(lamsin1(10) ,10) - 0.5 , 0.5 -0.5 ], [OUTatanvfinal(lamsin1(10) ,10) + 0.5 , 0.5 +0.5 ],
color=cap, label=false , lw=lwa )
Plots.plot!( [OUTatanvfinal(lamsin2(10) ,10) - 1.5 , 0.5 -1.5 ], [OUTatanvfinal(lamsin2(10) ,10) + 1.5 , 0.5 +1.5 ], color=cap,
label=false, lw=lwa )
Plots.plot!([OUTatanvfinal(lamsin2(10) ,10) - 1.5 , linapp_v2(9.999) - linapp_u2(9.999) ],
[OUTatanvfinal(lamsin2(10) ,10) + 1.5 , linapp_v2(9.999) + linapp_u2(9.999) ], color=cap, label=false, lw=lwa)
Plots.plot!( [OUTatanvfinal(lamsin1(10) ,10) - 0.5 , linapp_v(9.999) - linapp_u(9.999) ],
[OUTatanvfinal(lamsin1(10) ,10) + 0.5 , linapp_v(9.999) + linapp_u(9.999) ], color=cap, label=false , lw=lwa )
#### APPARENT GORNY
cap2=:cyan3
Plots.plot!( x0-> linupav1(x0) -linupau1(x0), x0-> linupav1(x0) +linupau1(x0), up_x0ran, color=cap2,
label=L"\mathrm{second} \, \mathrm{AH}", lw=lwa )
Plots.plot!( x0-> linupavB(x0) -linupauB(x0), x0-> linupavB(x0) +linupauB(x0), up_x0ran2, color=cap2, label=false, lw=lwa )
Plots.plot!( x0-> linupavBW(x0) -linupauBW(x0), x0-> linupavBW(x0) +linupauBW(x0), up_x0ran3, color=cap2, label=false, lw=lwa )
Plots.plot!( [linupav1(st)-linupau1(st) , linupavB(st)-linupauB(st) ], [linupav1(st)+linupau1(st) ,linupavB(st)+linupauB(st) ],
color=cap2, label=false, lw=lwa )
Plots.plot!( [linupavBW(przec2)-linupauBW(przec2) , linupavB(przec2)-linupauB(przec2) ],
[linupavBW(przec2)+linupauBW(przec2) ,linupavB(przec2)+linupauB(przec2) ],
color=cap2, label=false, lw=lwa )
Plots.plot!( [0.5-OUTatanufinal3(lamsin3(10),10) , 0.5-1.5 ],[0.5+OUTatanufinal3(lamsin3(10),10) , 0.5+1.5 ],
color=cap2, label=false, lw=lwa )
Plots.plot!( [1.5-OUTatanufinal3(lamsin4(10),10) , 1.5-1.5 ],[1.5+OUTatanufinal3(lamsin4(10),10) , 1.5+1.5 ],color=cap2,
label=false, lw=lwa )
Plots.plot!( [0.5-OUTatanufinal3(lamsin3(10),10) , linupav1(9.999)-linupau1(9.999) ],
[0.5+OUTatanufinal3(lamsin3(10),10) , linupav1(9.999)+linupau1(9.999) ], color=cap2, label=false, lw=lwa )
Plots.plot!( [1.5-OUTatanufinal3(lamsin4(10),10) , linupavBW(9.999)-linupauBW(9.999) ],
[1.5+OUTatanufinal3(lamsin4(10),10) , linupavBW(9.999)+linupauBW(9.999) ],color=cap2, label=false, lw=lwa )
### HORYZONTY ZEWNETRZNE
Plots.plot!( [0.5-0.5, 0.5-1.5 ], [0.5+0.5, 0.5+1.5 ], color=:grey, label=L"x_- \, \mathrm{horizons}" )
Plots.plot!( [0.5-1.5, 1.5-1.5 ] , [0.5+1.5, 1.5+1.5 ], color=:grey, label=false )
######### NULL INFINITIES
nl=:darkgreen
Plots.plot!( [-0.5+0.5, 0.5+0.5], [-0.5-0.5, 0.5-0.5], color=nl, label=L"\mathrm{null} \, \mathrm{infinities}" )
Plots.plot!( [0.5-(-0.5) ,0.5-(0.5) ], [0.5+(-0.5) ,0.5+(0.5) ], color=nl, label=false )
Plots.plot!( [1.5-1.5, 2.5-1.5], [1.5+1.5, 2.5+1.5], color=nl, label=false )
Plots.plot!( [ 2.5-1.5, 2.5-2.5 ], [ 2.5+1.5, 2.5+2.5 ], color=nl, label=false )
### TIMELIKE SINGULARITY
x1=0
Plots.plot!( t-> vsin(t,x1)-usin(t,x1), t-> vsin(t,x1)+usin(t,x1), range(-15000,15000,length=4), color=:red2, ls=:dot,linewidth=1,
label=L"{\mathrm{singularity}}")
######## POW POCZATKOWA
Plots.plot!( x0-> innv(x0)-innu(x0), x0-> innv(x0)+innu(x0), intabx0, color=:black,
label=L"\mathrm{initial} \, \mathrm{surface}", ls=:dashdot )
######## TRAJEKTORIE
lw1=0.6
Plots.plot!( l -> vtrajB(l) - utrajB(l),l -> vtrajB(l) + utrajB(l),tablamBb , aspectratio=1, color=:blue, label=L"x_0=5", lw=lw1 )
Plots.plot!( l -> vtrajC(l) - utrajC(l),l -> vtrajC(l) + utrajC(l),tablamCc , aspectratio=1, color=:green, label=L"x_0=0", lw=lw1 )
Plots.plot!( l -> vtrajD(l) - utrajD(l),l -> vtrajD(l) + utrajD(l),tablamDd , aspectratio=1, color=:orange,
linewidth=1, linestyle=:dash, label=L"x_0=11", lw=lw1)
Plots.plot!( l -> vtrajE(l) - utrajE(l),l -> vtrajE(l) + utrajE(l),tablamEe , aspectratio=1, color=:brown,
linewidth=1, linestyle=:dash, label=L"x_0=15", lw=lw1)
Plots.plot!( l -> vtrajA(l) - utrajA(l),l -> vtrajA(l) + utrajA(l), tablamAa, aspectratio=1, color=:red, label=L"x_0=10", lw=lw1 )
Out[37]:
In [38]:
###### ATRYBUTY
Plots.plot!(xlabel=L"\tilde{v} - \tilde{u}",
ylabel=L"\tilde{v} + \tilde{u}")
Plots.plot!(xlims=(-4.25,3.75), ylims=(-1.25,5.25), grid=false, xticks=-3:4:1, yticks=-1:6:5, legendfontsize=6,
xtickfontsize=6,ytickfontsize=6, xguidefontsize=10,yguidefontsize=10 )
######### LENS
lens!([-2.2, -1.7], [1.6, 2.1]; inset=(1, bbox(0.1, 0.2,0.21, 0.21)),
subplot=2, ticks=false, framestyle=:box, ls=:dash, lw=0.5 , lc=:grey, aspectratio=1)
lens!([-1.1,-0.6 ], [-0.1, 0.4]; inset=(1, bbox(0.15, 0.72,0.21, 0.21)),
subplot=2, ticks=false, framestyle=:box, ls=:dash, lw=0.5 , lc=:grey, aspectratio=1)
Out[38]:
In [39]:
Plots.savefig( "dagINNY.svg")
Out[39]:
"C:\\Users\\mihalbobula\\Desktop\\ew_final_przec\\dagINNY.svg"
In [ ]: